This page was rendered on: 2023-04-14 17:11:39

Introduction

This RMarkdown file perform analyses on the meteorological data that will be used in the DART (Dengue Advanced Readiness Tools) project. The meteorological data for this project is provided by ECMWF and USTH.

About the data

ECMWF

• Temporal coverage: 02/01/2019 - 31/12/2019

• Temporal resolution:

  • 6-hour intervals
  • generated every Monday and Thursday
  • 10 days forecasts from each start date

• Spatial coverage: Vietnam (100-110°E, 5-25°N)

• Spatial resolution: 0.25°

• 11 ensemble members

• Variables (unit):

  • t2m - (K) Temperature at 2m above the surface
  • tp - (m) total precipitation within the last 6 hours (e.g. timestamp of 12:00 means total precipitation between 0600 and 1200 – note this means that value at t=0 is always 0)

• Description: …

USTH

• Temporal coverage: 02/08/2018 - 30/08/2018

• Temporal resolution: hourly

• Spatial coverage:

  • Ho Chi Minh City: 105-108°E, 9-12°N

  • Hanoi: 104-107°E, 19-22°N

• Spatial resolution: 2km

• Variables (unit):

  • T2m - (°C) Temperature at 2m above the surface
  • Humidity - (%) Relative humidity at 2m above the surface

• Description: Downsampled and interpolated with WRF.

Library imports

library(tidyverse)
library(lubridate)
library(magrittr)
library(ncdf4)
library(sf)
library(stars)
library(gganimate)

Geospatial data

vn_lvl2_shp <- st_read("geospatial/gadm41_VNM_2.shp")

## Reading layer `gadm41_VNM_2' from data source 
##   `C:\Users\tuyenhn\Documents\GitHub\Meteorology\geospatial\gadm41_VNM_2.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 710 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 102.1446 ymin: 8.381355 xmax: 109.4692 ymax: 23.39269
## Geodetic CRS:  WGS 84

Read NetCDF files

• Here we will create a function that can read and generalized the NetCDF files from both ECMWF and USTH. The function will read the files, extract out all the variables and dimensions into a named list.

ingest_ncdf <- function(fpath) {
  # read NetCDF file and get list of vars and dims
  ncdf4_temp <- nc_open(fpath)
  varlist <- as.list(c(names(ncdf4_temp$var), names(ncdf4_temp$dim)))

  # extract the vars and dims into a named list
  ncdf_obj <- varlist %>%
    map(\(x)
    if (x %in% names(ncdf4_temp$var) ||
      (x %in% names(ncdf4_temp$dim) && ncdf4_temp$dim[[x]]$dimvarid$id != -1)) {
      ncvar_get(ncdf4_temp, x, verbose = FALSE)
    })
  names(ncdf_obj) <- varlist

  nc_close(ncdf4_temp)

  # remove NULL vars and dims
  ncdf_obj %>% compact()
}

# example
usth_t2_hcm <- ingest_ncdf("USTH-WRF/T2_HCM.nc")

ecmwf_01_02 <- ingest_ncdf("ECMWF/01-02-surface.nc")

• With the data extracted, we can see what are the variables

# extracted var names
names(usth_t2_hcm)

## [1] "latitude"  "longitude" "time_data" "T2m"       "Humidity"

names(ecmwf_01_02)

## [1] "t2m"    "tp"     "time"   "lon"    "lat"    "number"

Transform to tibble dataframe

• Here we will create another function that will transform this named list into a tibble (dataframe). Each row of this tibble will represent the weather variable in a specific point in time at a specific longitude and latitude. This function assumes that the inputted object is a named list as exported from ingest_ncdf(); or of similar structure.
• Why? By transforming into a tibble, the data can then be used to plot figures and animations with ggplot2. The tibble will only have 4 columns: Latitude, Longitude, Weather variable, and Time variable. This also requires the inputted object to have exactly 4 of those dimensions as well (i.e. you should aggregate ensemble members into a single layer).
• Function specs: TBC

Please note that this is a memory intensive process. Depends on the total amount of pixels, the final tibble size can be upwards of 400MB!

ncdf_to_tibble <- function(lon, lat, var, time, time_start, time_step) {
  res_df <- tibble(lon = double(), lat = double(), var = double())

  time_ticks <- length(time)

  pb <- progress::progress_bar$new(total = time_ticks)
  for (i in 1:time_ticks) {
    pb$tick()
    temp_res_df <- tibble(
      lon = as.vector(lon),
      lat = as.vector(lat),
      var = as.vector(var[, , i]),
      time = ymd_hms(time_start, tz = "Asia/Ho_Chi_Minh") +
        time_step(time[i])
    )
    res_df %<>% bind_rows(temp_res_df)
  }

  res_df
}

# example
# usth_t2_df <- ncdf_to_tibble(
#   lon = usth_t2_hcm$longitude,
#   lat = usth_t2_hcm$latitude,
#   var = usth_t2_hcm$T2m,
#   time = usth_t2_hcm$time,
#   time_start = "2018-08-02 00:00:00",
#   time_step = lubridate::hours
# )

# ECMWF's data requires a bit more transformation to fit the function
# ecmwf_df <- ncdf_to_tibble(
#   lon = replicate(length(ecmwf_01_02$lat), ecmwf_01_02$lon),
#   lat = t(replicate(length(ecmwf_01_02$lon), ecmwf_01_02$lat)),
#   var = apply(ecmwf_01_02$t2m, c(2, 3, 4), mean), # aggregate ensemble members by mean
#   time = ecmwf_01_02$time,
#   time_start = "2019-01-02 00:00:00",
#   time_step = lubridate::hours
# )

# it is recommended to only run the function once then save it as a tibble for later use
# write_rds(usth_t2_df, "USTH_T2_HCMC.rds")
# write_rds(ecmwf_df, "ECMWF_2019_01_02.rds")
usth_t2_df <- read_rds("USTH_T2_HCMC.rds")
ecmwf_df <- read_rds("ECMWF_2019_01_02.rds")

Visualization

• We can now plot the raster at one point in time or produce an animation. Since the only change between figures is the data itself, we will create a function with predefined styling and elements.
• Note that animations rendering might take longer than expected.

Plotting function

plot_weather <- function(df) {
  ggplot(data = df) +
    geom_sf(data = vn_lvl2_shp) +
    geom_raster(aes(x = lon, y = lat, fill = var), alpha = 0.8) +
    coord_sf(
      xlim = c(min(df$lon), max(df$lon)),
      ylim = c(min(df$lat), max(df$lat))
    ) +
    scale_fill_gradient(low = "green", high = "red")
}

ECMWF t2m static

# filter to get only the first time step
ecmwf_t2_static <- ecmwf_df %>%
  filter(
    time == ecmwf_df %>%
      pull(time) %>%
      head(1)
  )

ecmwf_t2_static %>%
  plot_weather()

ECMWF animation

gganim_p <- plot_weather(ecmwf_df) +
  transition_time(time) +
  labs(title = "{frame_time}")

animate(gganim_p, fps = 5)

anim_save("t2m_vietnam_ECMWF_2019-01_5fps.gif")

USTH t2m static

# filter to get only the first time step
usth_t2_static <- usth_t2_df %>%
  filter(
    time == usth_t2_df %>%
      pull(time) %>%
      head(1)
  )

usth_t2_static %>%
  plot_weather()

USTH animation

gganim_p <- plot_weather(usth_t2_df) +
  transition_time(time) +
  labs(title = "{frame_time}")

animate(gganim_p, fps = 5)

anim_save("t2m_hcmc_USTH_2018_5fps.gif")